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We use our newly constructed quantum Gibbs ensemble Monte Carlo algorithm to perform com¬ 
puter experiments for the two phase coexistence of a hydrogen-helium mixture. Our results are in 
quantitative agreement with the experimental results of C. M. Sneed, W. B. Streett, R. E. Sonntag, 
and G. J. Van Wylen. The difference between our results and the experimental ones is in all cases 
less than 15% relative to the experiment, reducing to less than 5% in the low helium concentration 
phase. At the gravitational inversion between the vapor and the liquid phase, at low tempera¬ 
tures and high pressures, the quantum effects become relevant. At extremely low temperature and 
pressure the hrst component to show superhuidity is the helium in the vapor phase. 
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I. INTRODUCTION 

Hydrogen and helium are the most abundant elements 
in the Universe. They are also the most simple. At am¬ 
bient conditions helium is an inert gas with a large band 
gap. Because of its low mass and weak inter-atomic in¬ 
teractions, it has fascinating properties at low tempera¬ 
tures such as superfluidity. The molecular hydrogen and 
helium mixture is therefore of special theoretical impor¬ 
tance since it is made by the two lightest elements in 
nature which have the lowest critical temperatures. This 
mixture is found to make the atmosphere of giant planets 
like the Jovian and is essential in stars. 

An important problem to study is the phase coexis¬ 
tence of the fluid mixture and the determination of its 
coexistence properties. Some early experimental studies 
n have shown that at coexistence, at low tempera¬ 
ture, the mixture shows a strong asymmetry in species 
concentrations in the liquid relative to the vapor phase, 
with an abundance of helium atoms in the vapor. This 
phenomenon results in the liquid floating above its vapor 
[3| since helium has approximately twice the molecular 
weight of hydrogen. Such experimental coexistence stud¬ 
ies has later been extended at higher temperature and 
pressure Si allowing to determine a quite complete 
picture for the coexistence phase diagram of this mixture 
in the temperature range from 15.5 K to 360 K and in 
the pressure range from 5 bars to 75 kbars. Another in¬ 
teresting issue is whether this system exhibits fluid-fluid 
solubility at extremely high pressure [§-[T^. a situation 
hard to achieve in the laboratory. 

In this work we perform a numerical experiment for 
the two phase coexistence problem of the hydrogen- 
helium mixture at low temperatures and pressures using 
the Quantum Gibbs Ensemble Monte Carlo (QGEMC) 
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method recently devised [H, [3 to solve the coexistence 
of a generic quantum boson fluid where the particles 
interact with a given effective pair-potential. We will 
be concerned with situations where the absolute tem¬ 
perature, T, and the number density, pa, of each one 
of the two components a = a, b oi mass TOq,, are such 
that at least one of the two components is close to its 
degeneracy temperature (Td)o, = /makB, with 

ks Boltzmann constant. Eor temperatures much higher 
than max{(rD)a} quantum statistic is not very impor¬ 
tant. This path integral Monte Carlo simulation enables 
us to study the quantum fluid mixture from first princi¬ 
ples, leaving the effective pair-potentials between the two 
species, the hydrogen molecules and the helium atoms, as 
the only source of external information. There are stud¬ 
ies on reproducmg such coexistence from an equation of 
state approach [l^. Our QGEMC method is expected to 
break down at high densities near the solid phase. More¬ 
over, clearly our approach becomes not anymore feasible 
at extremely high pressures when the hydrogen is ionized 
and one is left with delocalized metallic electrons nn. 

Our binary mixture of particles, of two species la¬ 
beled by a Greek index, with coordinates R = = 

1 , 2 ,..., Na and a = a, b}, and interacting with a central 
effective pair-potential (j)a/ 3 {f), has a Hamiltonian 

2 Af„ 

Ct — l ia — 'i 

^ ( 1 ) 

a,P=l icjfi 

where the prime on the sum symbol indicates that we 
must exclude the terms with ia = jp when a = j3 and 
Act = k? 12m a- 

The density matrix for the binary mixture at equilib¬ 
rium at an absolute temperature T is then p = with 

/3 = l/fcsT. Its coordinate representation p{R,R',j3) 
can be expressed as a path (i?(r)) integral in imagi- 



2 


nary time (r) extending from R = i?(0) to R' = R{fi) 
[a. The many-particle path is made of = Na + Ni, 
single-particle world-lines which constitute the configu¬ 
ration space one needs to sample. Since the Hamilto¬ 
nian is symmetric under exchange of like particles we can 
project over the bosonic states by taking pb{R, R', 13) = 
p[R,VR',l3)/{Na\Nb\) where V indicates a permu¬ 
tation of particles of the same species. 

If we call p the number density of the mixture, Xa the 
molar concentration of species a (xb = 1 — Xa), P = 
P{T,p,Xa) the mixture pressure, and Pa = l^a{T,P,Xa) 
the chemical potential of species a, we want to solve the 
two phase, I and II, coexistence problem 

PaiT,P,xi^^)=Pa{T,P,xi^^^) (2) 

PbiT,P,xi^^)=Pb{T,P,xi^^^) (3) 


for the concentrations, xi^'^ and x^^\ (and the densities, 
pP') and pPP) in the two phases. Since our mixture is 
not symmetric under exchange of the two species, a and 
b, we expect in general 1 — x[!^\ 

Our QGEMC algorithm [ij uses two boxes maintained 
in thermal equilibrium at a temperature T and contain¬ 
ing the two different phases. It employs a menu of seven 
different Monte Carlo moves: the volume move (g = 1) 
allows changes in the volumes of the two boxes assuring 
the equality of the pressures between the two phases, the 
open-insert (g = 2), close-remove (g = 3), and advance- 
recede (g = 4) allow the swap of a single-particle world¬ 
line between the two boxes assuring the equality of the 
chemical potentials between the two phases, the swap 
(g = 5) allows to sample the particles permutations, 
and the wiggle (g = 6) and displace (g = 7) to sam¬ 
ple the configuration space. We thus have a menu of 
seven, g = 1, 2,..., 7, different Monte Carlo moves with 
a single random attempt of any one of them occurring 
with probability Gg = g,/ V =i 9g- 

The paper is organized as follows: In Section|n]we de¬ 
scribe the particular binary mixture studied; in Section 
m we describe the simulation method employed; in Sec¬ 
tion |IV] we present our numerical results; Section |V] is for 
final remarks. 


a hydrogen molecule and a helium atom is the Roberts 
[Sill. All can be put in the following central form 


^(r) = £<i)(x) 
4)(a;) = exp(a 



— fix — "fx^) — 


Cs 


Cm\ 

) 


F{x), 


f exp[—(H/a: — 1)^] x<D 
1 X > D ^ 


( 6 ) 

(7) 

( 8 ) 


where x = r/rm, with the position of the mini¬ 
mum, and the various parameters are given in Table HI 
We have </'HeHe(’’"i) = -10-8, = -34.3, and 

<^112 He (’""i) = —14.8. Moreover we have a slight pos¬ 
itive non-additivity: [rf^]H 2 Ue = 3.375 > ([rj^JneHe + 
[’’m]H2H2)/2 = 3.189. 

The experimental coexistence data [I|, [p, [3 is given 
in Table I of the supplemental material [2ll| and rep¬ 
resented schematically in Fig. [TJ For example, the 



FIG. 1. (color online) Schematic pressure-composition phase 
diagram for six isotherms of the hydrogen-helium mixture at 
low temperatures and low pressures as drawn in the experi¬ 
mental work by Streett at al. 0]. The thick continuous black 
line is the mixture critical line. 


II. THE Ha-HE MIXTURE 

We consider a binary fluid mixture of molecular hydro¬ 
gen (H 2 ) and the isotope helium four (^He), two bosons. 
We take 1 A as unit of lengths and ks K as unit of en¬ 
ergies. We indicate with an asterisk over a quantity its 
reduced adimensional value. We have for the parameter 
Act = bi/2ma of the two species a = H 2 , ^He 

= 12.032, (4) 

Ane = 6.0596. (5) 

The pair-potential between two helium atoms is the Aziz 
et al. [13 HFDHF2, the one between two hydrogen 
molecules is the Silvera et al. [I^, and the one between 


mixture at T = 31 K has a lower critical state at 
P = 0.207 /cbKA” , a;He = 0.214 and an upper criti¬ 
cal state at P = 1.96 /cbKA ^,a:He = 0.49. The set of all 
critical states constitutes the a;—line, T = T^fP), such 
that for T > Tx then = x^l\ The experimental 

X-line of Sneed et al. is shown in Fig. [5] for the 

low temperature and low pressure mixture. In the figure 
we also show the experimental line for the gravitational 
inversion described in Section IIII Al 

For temperatures higher than the hydrogen critical 
point rH 2 = 33.19 K (Phs = 0.094 fc^KA there is 
only an upper critical point [^. On the temperature at 
which Tx{P) reaches its minimum there is no unanimous 
consensus among the various experimental works. 
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TABLE I. Pair-potentials parameters: (/>pair 


pair 


'^m 

a 

P 

7 

Ge 

Cs 

Gio 

D 

He-He 

10.8 

2.9673 

13.208 

13.353 

0 

1.3732 

0.42538 

0.17810 

1.2413 

H 2 -H 2 

315778 

3.41 

1.713 

10.098 

0.41234 

1.6955x10"'^ 

7.2379x10"® 

3.8984x10"® 

1.28 

H 2 -He 

14.76 

3.375 

13.035 

13.22 

0 

1.8310 

0 

0 

0.79802 



PlksKA-®) 


FIG. 2. (color online) Reproduction of Fig. 3 of Sneed et al. 
0 for the X-line and the g-line (see Eq. fTTl))- The inset 
shows the two lines in the temperature-composition plane. 


III. SIMULATION METHOD 


We use our QGEMC method, described in Ref. 
where we monitor the number densities of the two co¬ 
existing phases, = (A^^] + 

with i = I, II, the concentrations of He in the two 
phases, = .^He/(^He + -^h]) < 1 ; the pressure 
P. We shall conventionally order so that / 

will be the vapor phase and II the liquid phase, unless 
pP) = pHt) ^ jjj which case we have a fluid-fluid phase co¬ 


existence. In the simulation we fix: N = 


{II) 


N, 


(I) 


with 

ti2 rL2 


N, 


{II) 


Ha “^"'Ha — XW} 


{!) 

He 


and 


V = -h . otherwise iV^e > ^ > ^4a^^ and 

are allowed to fluctuate keeping 
and -|- constants. The Gibbs 

phase rule for a two phase coexistence of a binary mixture 
assures that one has two independent thermodynamic 
quantities 22|. So our control parameters will be the 
absolute temperature T and the global number density 
p = N/V (instead of the pressure as in the experimental 
case). As usual a finite N sets the size error for our calcu¬ 
lation. Whereas x > 0 will regulate the size asymmetry 


numerical effect so that for 


-x^g^(l-f X)] 




{II), 


-'''He ~ 2 ^ 


(1 + x)(a'Hi - 4^/^) 


> 0 , 


N 


■X 


-<e > 0 , 


(9) 

( 10 ) 


if must have 0 < 4g^i < 1/(1 -I- x) < 

4g < 1 and if < x^g^^, then 0 < x^g < 1/(1 + 
x) < a:^gi < 1. Moreover we must also always have 
< p < p^^^\ The initial condition we chose for our 
simulations was always as follows: p^-^i = = p and 

4i = 44 = 1/(1 + X)- 

Due to the short-range nature of the effective pair- 
potentials of Eq. (|ni) we will approximate, during the 
simulation, ^(r) = 0 for r > [j’mjHaHa (this 

corres pon ds to the truncated and not shifted choice in 
Ref. |23l|). Where in order to comply with the mini¬ 
mum image convention for the potential energy calcula¬ 
tion, we make sure that the conditions > 2rcut, 

for i = I, II, are always satisfied during the simula¬ 
tion. This approximation is the only other source of 
error apart from the size one. The two are related be¬ 
cause for instance in the fluid-fluid coexistence, when 
I/(^) ?a l/G-f) ~ yj2 during the simulation, we require 
r-cut ~ {Nl2pY/^l2 » [r^JusHs for some given p. 

The path integral discretization imaginary time step 
5t = P/K, with K the number of time slices, is cho¬ 
sen so that 5t* = 0.002, which is considered sufficiently 
small to justify the use of the primitive approximation 
of the inter-action [T^. The parameters M, defined in 
[T3| . will be called Mq for each relevant move q and 
the parameter An, also defined in [T^, is always cho¬ 
sen equal to 0.01. In order to fulfill detailed balance we 
must choose M 2 = Ms- In particular we always chose 
M 2 = 5, M 3 = 5, M 4 = 5, Ms = 5,M^6 = 5. Regarding 
the frequency of each move attempts, we always chose 
gi = 0.001,32 = 1,33 = 1,54 = 1,35 = l,ff6 = 1,37 = 
0.1. The parameter C defining the relative weight of the 
Z and G sectors [3 is adjusted, through short test runs, 
so as to have a Z sector frequency as close as possible 
to 50%. We accumulate averages over 10® blocks each 
made of 10 ® attempted moves with quantities measured 
every 10® attempts. Since the volume move is the most 
computationally expensive one we chose its frequency as 
the lowest. During the simulation we monitor the accep¬ 
tance ratios of each move. The various simulations took 
no more than ~ 150 GPU hours on a 3 GHz processor. 
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A. Barotropic phenomenon and gravitational 
inversion 

The condition for the gravitational inversion observed 
experimentally Q is 

> 

, (H) 

where mHe/™H 2 = 1-98553. When this condition on the 
mass density inversion respect to the number density is 
satisfied, the liquid phase will float on top of the vapor 
phase. The condition of Eq. (HU can also be rewritten 
as 

(l + , ( 12 ) 

where k = raYitilrriYi^ — 1 = 0.98553. This condition may 
be satisfied when the concentration of He in the vapor 
phase is bigger than in the liquid phase, at low tempera¬ 
tures, and the number density of the liquid is close to the 
one of the vapor, at high pressures. We expect quantum 
effects to become important in this regime, before solid¬ 
ification which is expected to occur for T < Ts{P). The 
gravitational inversion of Eq. (HU will be satisfied for 
T < Tg{P). The experimental s-line T = Ts{P) and g- 
line T = Tg{P) have been determined in Fig. 3 of Sneed 
et al. 0 ) in the laboratory. 


B. Pressure calculation 

We will use the virial estimator for the pressure (see 
Eq. (6.18) of Ref. 0 ). With long-range corrections [2^ 
which can be quite big in the liquid phase. More details 
on the pressure calculation are given in the supplemental 
material [ 2 l|. 


IV. NUMERICAL RESULTS 

Our results are summarized in Table HIl and compared, 
in Fig. [3 with the experimental data of Refs. 0,1 
Wl summarized in a Table in the supplemental material 

13). 

In all studied cases we chose N = 128 and St* = 0.002. 
We explored the vapor-liquid coexistence (in this work 
we will denote with “vapor-liquid” coexistence one where 
^ five temperatures, T = 2,5,15.5,26,31 

degrees Kelvin, and the fluid-fluid coexistence (in this 
work we will denote with “fluid-fluid” coexistence one 
where at T = 31 K. For the first two lower 

temperatures studied we could not find any experimental 
data for a comparison. In these two cases when we put a 
number with trailing dots in the table it means that after 
the initial equilibration period the measured property did 
not change anymore during the rest of the simulation. 



^He 


FIG. 3. (color online) Comparison between the results of our 
numerical experiments, points from Table El and of the lab¬ 
oratory experiments, lines from Table I in the supplemental 
material [^, for the pressure-composition of three isotherms 
of the hydrogen-helium mixture phase diagram. A logarith¬ 
mic scale is conveniently used on the ordinates. The double 
circled points at T = 15.5 K denote the case where we ob¬ 
served gravity inversion in the numerical experiment. 


For the temperature T = 15.5 K, as it can be readily 
verified using the relation of Eq. HU, we observe gravi¬ 
tational inversion on the point at p = 0.02 A when the 
component with the highest degeneracy temperature is 
the hydrogen in the liquid phase with Tjj Ri 2 K. Clearly 
choosing higher pressures quantum statistics will become 
more and more important for the fluid mixture before 
reaching the solid state. 

For the points at T = 26 K, T = 31 K,p = 
0.006 A“^ X = 116/12, and T = 31 K, p = 0.03 A“^ y = 
1 we observed exchanges of identity between the two 
phases, during the simulation. 

At a temperature T = 31 K and a pressure of P = 
0.07(2) keKk we found a vapor-liquid coexistence, 
choosing y = 116/12. This point should be subject to 
greater size error than all other points simulated, and 
be thus the less reliable, since we only have, in the two 
boxes, a total of 12 helium atoms. Increasing the pres¬ 
sure to P = 0 . 21 ( 2 ) ^bKA in agreement with the 
experiment, we did not find coexistence and we observed 
p(^) Ri pPP PS p and x^l ~ x^^} ~ 1/(1 -by). Increasing 

the pressure to P = 3.5(4) A:bKA , we did not observe 
exactly p*-^^ = p^^^\ as measured in the fluid-fluid tran¬ 
sition observed in the laboratory 0|. The same holds 
true for the point at the same temperature but higher 
pressure P = 6.3(6) fc^KA 

For all measured points except the one at the lower 
temperature, P = 2 K of Table UTl the superfluid fraction 
of the two components in either phase was negligi¬ 
bly small. At T = 2 K of Table [TTl below the helium 
lambda-temperature, we observed a negligible superfluid 
fraction of both components in the liquid phase and of 
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the hydrogen in the vapor phase. The helium in the vapor 
phase was found to have a superfluid fraction of 0.012(3), 
indicating a tendency to supefluidity. 

When we do not observe exchanges of identity between 
the two phases, during the simulation, we are able to 
find accurate average values for the various measured 
quantities. Otherwise a histogram analysis of the data 
is necessary with a non-linear fit using the superposi¬ 
tions of two shifted Gaussians. For example in Fig. 0] 
we show the procedure used to extract the helium con¬ 
centrations of the two coexisting phases for the case 
T = 26 K, p = 0.01 X = 90/38. 

4.0 
3.5 
3.0 

a 2.5 

cc 

^ 2.0 
-s 

'■S 1.5 

1.0 
0.5 
0.0 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 


FIG. 4. Fit of the histogram for the block averages of both 
with the sum of two gaussians with six param¬ 
eters. This is case T = 26 K,p = 0.01 A ^,x = 90/38, where 
we had box identity excanges. 

The measured property which is less accurate is the 
pressure due to the size error and the long-range correc¬ 
tion dependent on the Tcut choice. This problem could be 
overcome by using instead of the N, V, T version of the 
Gibbs ensemble algorithm its N, P, T one [2^ . 

A. Finite size effects 

We studied the finite size effects atr = 31K,p = 
0.03 A ^,x = 1. In Table IIIII we show the results 
for the isothermal pressure-composition coexistence at 
N = 64,128, and 256. As the number of particles in¬ 
creases we observe a decrease in the ratio of number of 
exchanges of identity between the two phases and total 
number of particles: For = 64 the exchanges occurred 
many times, for N = 128 only once, and for N = 256 
never. For the case A = 64 we found the peak of the 
first gaussian for the hystogram of xue with a negative 
value. The simulation with A = 64 took 1.0 x lO^s, the 
one with A = 128 took 1.6 x lO^s, and the one with 
A = 256 took 4.0 x 10®s. From the comparison we see 
how there is not much difference between A = 128 and 
A = 256. Apart from the smaller statistical errors in the 


latter case, the concentrations slightly differ in the two 
cases. 


B. Importance of the particle exchanges and of the 
quantum effects 

Setting to zero the frequency of the swap move at¬ 
tempts our algorithm reduces to a path integral calcu¬ 
lation for distinguishable particles obeying to the Boltz¬ 
mann statistics. On the other hand, choosing K = 2 
(with Mq = 1 for all q) and = AJjg 0 it reduces 
to the classical Gibbs Ensemble Monte Garlo (GEMC) 
algorithm of Panagiotopoulos . 

Eor the state point T = 15.5 K and p = 0.02 A ^ with 
A = 128, we performed two simulations for each of the 
two cases suggested above to estimate the importance 
of particles exchanges which underlies the Bose-Einstein 
statistics and of quantum effects, respectively. To reach 
the GEMG limit from our QGEMG algorithm we chose, 
in particular, A^^ = AJjg = 10“^. The results are shown 
in Table HVl The acceptance ratio for the swap move was 
around 0.5 in the full quantum case and imposed zero in 
the other two simulations. 

As we can see from the table, for this state point, 
there is a very small difference between the path inte¬ 
gral simulation with the full Bose-Einstein statistics and 
the one with the Boltzmann statistics. In particular, only 
the densities of the vapor phase are different in the two 
cases. In both these simulations we observe the gravita¬ 
tional inversion. We expect that increasing the pressure 
and thereby the density or reducing the temperature the 
particles exchanges will become increasingly important. 

On the other hand, there is a large difference between 
these two simulations and the classical GEMG one. In 
particular, the gravitational inversion is not observed in 
the classical limit simulation, even if after a short equili¬ 
bration time the simulation converged towards the condi¬ 
tion = 1, i.e. all helium atoms, the heaviest species 
in the mixture, were found in the less dense phase. 

V. CONCLUSIONS 

In conclusion, we performed path integral Monte Carlo 
simulations, using our newly developed QGEMG method, 
for the two phase coexistence of the hydrogen-helium 
mixture away from freezing. This asymmetric mixture 
displays at low temperature, a big concentration asym¬ 
metry in the two coexisting phases, whereas the densities 
of the two phases tend to become equal at high pressure. 
This is responsible for a gravitational inversion, where 
the liquid, the more dense phase, with an abundance of 
hydrogen, floats above the vapor, the less dense phase, 
with an abundance of helium. In this coexistence region 
of the temperature-pressure diagram, quantum statistics 
is expected to play an important role and in our simula¬ 
tions we are able to observe such gravitational inversion. 
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TABLE II. Numerical isothermal pressure-composition at coexistence. We always used N = 128 and St* = 0.002. 


T(K) 

P (A-^) 

X 

P (fcsKA"®) 

*^He 


(A-^) 

pW (A-^) 

2.0 

0.015 

1 

-0.08(7) 

0.214. .. 

0.639... 

0.02456(1) 

0.012605(2) 

5.0 

0.010 

1 

0.014(2) 

0.1787(1) 

1.000... 

0.025910(6) 

0.005113(1) 

15.5 

0.010 

1 

0.093(7) 

0.00457(9) 

0.948(1) 

0.02410(1) 

0.006544(5) 

15.5 

0.015 

1 

0.16(4) 

0.0125(3) 

0.923(1) 

0.02304(2) 

0.011525(7) 

15.5 

0.020 

1 

0.30(9) 

0.0142(4) 

0.921(1) 

0.02373(2) 

0.017619(5) 

26.0 

0.010 

90/38 

0.14(2) 

0.044(2) 

0.546(4) 

0.01890(5) 

0.00669(1) 

26.0 

0.015 

90/38 

0.25(5) 

0.118(3) 

0.593(8) 

0.01888(7) 

0.01105(5) 

26.0 

0.020 

90/38 

0.6(1) 

0.170(3) 

0.69(3) 

0.02115(2) 

0.01759(8) 

31.0 

0.006 

116/12 

0.11(3) 

0.091(1) 

0.201(7) 

0.014(2) 

0.00564(6) 

31.0 

0.008 

1 

0.21(2) 

0.5025(6) 

0.511(1) 

0.008016(7) 

0.00795(1) 

31.0 

0.030 

1 

3.5(4) 

0.832(4) 

0.113(3) 

0.03198(4) 

0.02805(1) 

31.0 

0.035 

1 

6.3(6) 

0.932(2) 

0.0243(9) 

0.03955(5) 

0.03111(1) 


TABLE III. Numerical isothermal pressure-composition coexistence at T = 31 K,p = 0.03 A = 1 as a function of the 
number of particles N. We always used St* = 0.002. 


N 

'f 'cut (-^) 

P (fcsKA ") 

(II) 

®He 


1 

to 

p(^) (A-^) 

64 

5 

2.4(8) 

0.83(3) 

- 

0.03144(7) 

0.02782(3) 

128 

6 

3.5(4) 

0.832(4) 

0.113(3) 

0.03198(4) 

0.02805(1) 

256 

8 

3.4(2) 

0.840(3) 

0.098(3) 

0.03180(3) 

0.028170(9) 


Our numerical experiments are also in good quantitative 
agreement with the experimental results of C. M. Sneed, 
W. B. Streett, R. E. Sonntag, and G. J. Van Wylen in the 
late 1960’s and early 1970’s. The difference between our 
results on the helium concentration in the two phases and 
the experimental ones is in all cases less than 15% in the 
high helium concntration phase and than 5% in the low 
helium concentration phase, relative to the experiment. 

These results for the hydrogen-helium mixture can be 
of interest for the study of cold exoplanets with an at¬ 
mosphere made predominantly by such a fluid mixture 
and with the right temperature and pressure conditions 
for there to be coexistence. In such cases it could be pos¬ 
sible to observe the gravitational inversion phenomenon 
and consequent changes in the planet moment of inertia, 
depending on the atmospherical and climatic conditions. 


At extremely low temperature and pressure we find 
that the first component to show superfluidity is the he¬ 
lium in the vapor phase. 

Our QGEMC method 0 is extremely simple to use, 
reduces to the Gibbs ensemble method of Panagiotopou- 
los [13 in the classical regime, and gives an exact numer¬ 
ical solution of the statistical physics phase coexistence 
problem for boson fluids. 

An open problem, currently under exam, is the influ¬ 
ence of the finite-size effects on the determination of the 
binodal curves close to the lower strongly asymmetric 
critical points, as for example in our case T = 31 K, p = 
0.006 A” ,X = 116/12. This requires additional simula¬ 
tions at an higher and lower number of particles. 
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